Long wavelength instability for uniform shear flow 
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Abstract 

Uniform shear flow is a prototype non-equilibrium state admitting detailed 

study at both the macroscopic and microscopic levels via theory and computer 

simulation. It is shown that the hydrodynamic equations for this state have 

a long wavelength instability. This result is obtained first from the Navier- 

Stokes equations and shown to apply at both low and high densities. Next, 

higher order rheological effects are included using a model kinetic theory. The 

results are compared favorably to those from Monte Carlo simulation. 
PACS numbers: 47.20.Ft, 05.20.Dd, 05.60.+W, 0.5.70Ln 
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In spite of extensive formal theory for non-equilibrium statistical mechanics, definitive 
tests and controlled illustrations outside the domain of linear response are rare. One such 
rare case is the macroscopic state of uniform shear flow. Like planar Couette flow the x- 
component of the average flow velocity varies linearly in the y direction. However, the flow 
is generated by periodic boundary conditions in the local Lagrangian frame leading 
to a uniform temperature and pressure, with a monotonic increase of the temperature. 
This is in contrast to planar Couette flow, driven by local boundary conditions, for which 
the temperature is non-uniform but stationary @]. The advantage of uniform shear flow 
is that the boundary conditions allow application of computer simulation methods with 
periodic imaging to emulate bulk effects for small systems, just as in the simulation of 
equilibrium states. In addition, there are many theoretical advantages such as the absence 
of a hydrodynamic boundary layer. The problems associated with increasing temperature 
can be controlled by the addition of a thermostat so that a steady state with the desired 
flow field results. During the past fifteen years uniform shear flow has been the focus of 
attention in most studies of classical non-equilibrium statistical mechanics by a wide range 
of theoretical and simulation methods It is one of the few means by which rheology 

and transport far from equilibrium can be studied at the fundamental level. More recently, 
simulations of uniform shear flow have provided the test for new concepts of non-equilibrium 
variational principles and for the relationship between transport and chaotic Hamiltonian 
flows g]. 

Although it has been known for more than ten years that uniform shear flow undergoes 
a transition to an ordered state at large shear rates 0, its stability at smaller shear rates 
has not been questioned. Our objective here is to report an instability of uniform shear flow 
at any shear rate for a sufficiently long wavelength perturbation. More precisely, solutions 
to the hydrodynamic equations linearized about the macroscopic state of uniform shear 
flow show exponential growth in time for wavenumbers smaller than a critical wavenumber 
kcr{ci) > for a > 0. For the case of zero shear rate these linear equations deflne the 
hydrodynamic modes which describe how small perturbations of the velocity, temperature, 
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and density fields decay back to equilibrium. These are the two sound modes, a heat mode, 
and two shear modes. In a similar way the hydrodynamic modes describing response to 
perturbations of any reference stationary state can be identified. If all such modes decay in 
time the reference stationary state is referred to as linearly stable; conversely, if one or more 
modes leads to an increasing amplitude in time the reference state is linearly unstable. 

The hydrodynamic equations are approximations to the more general local conservation 
laws for the mass, energy, and momentum densities, 

Dtn + nV -V^ 0, (1) 
Ae + (e + p) V • U + V • S + PijdiUj = w, (2) 

DtUi + p-^dip + p-^djP,j = 0, (3) 

where Dt = 9t+U- V is the material derivative. The momentum density is related to the flow 
field U by g = pU, p = mn is the mass density (m = mass, n = number density), and e is the 
internal energy density. The pressure, p, is not an independent variable but is defined to have 
the same functional relationship to n and e as in equilibrium. The inhomogeneous term on 
the right side of the energy equation is due to an external non-conservative force representing 
the thermostat. There are several thermostats that have been used in simulations and 
theory. Here the thermostat is fixed by a force on each particle proportional to the velocity 
relative to the local flow fleld. Consequently, its average value is zero and the effect of this 
force occurs only in the energy equation. The proportionality constant is determined by 
requiring stationarity for the uniform shear flow state. The resulting term w in Eq. (2) 
is the same as that obtained using the thermostat for simulations. Finally, the irreversible 
heat and momentum fluxes are denoted by S and Pij, respectively. These equations are 
exact but incomplete until the irreversible fluxes are specifled in terms of the hydrodynamic 
fields. Nevertheless, it is possible to define the state of uniform shear fiow without their 
detailed form. Steady state uniform shear flow is deflned by constant (in both space and 
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time) energy and mass densities, {eo,po}, and a flow velocity whose only non-vanishing 
component is Uq^ = ay. The constant, a, is the shear rate and provides the single control 
parameter measuring the deviation from equilibrium. The boundary conditions are simple 
periodic conditions in the local Lagrangian coordinate frame, r' = r — Uo(r)t. Substitution 
of these assumptions for e, p, and U into the above conservation laws shows they are all 
satisfied if S and Pij are also uniform and the parameter of the thermostat in the energy 
equation is chosen such that wq = aPoxy 

Consider small deviations of the hydrodynamic variables from the uniform shear flow 
state. To be explicit it is necessary to specify the heat and momentum fluxes. We first 
consider the case for which all spatial gradients are small, including the shear rate a. Then 
the heat flux is given by Fourier's law and the momentum flux is given by Newton's viscosity 
law, 

S = -AVT, (4) 

Pij = -vidiU, + d,U, - |V • V6,,) - kV ■ U5,„ (5) 

where A, t], and k, are the thermal conductivity, shear viscosity, and bulk viscosity, respec- 
tively. These are the leading terms in a "uniformity" expansion ordered according to spatial 
gradients of the hydrodynamic fields [0. It is convenient to choose density and temperature 
as independent thermodynamic parameters and to denote the deviations in the hydrody- 
namic parameters from their values for uniform shear flow by = {^n, 6T, 6Ui}. Also, we 
consider a fluid of hard spheres to specify the equations of state, e = e(n, T) and p = p{n, T). 



These are determined from the Percus-Yevick equation |Tl| which is known to be accurate 
over the entire fluid phase. Combining Eqs. (Ill)-® and retaining terms only up through 
linear order in identifies the hydrodynamic equations for small perturbations about uni- 
form shear flow. The boundary conditions are made explicit by looking for solutions of the 
form 

l/„(r, i) = ya(k, t) exp(ik ■ r'), (6) 



where ?/a(k, t) is the amphtude of a mode with wavelength, 211 /k. To simphfy the analysis 
and to focus on the instability the following is restricted to the special case — kr^ — 0, i.e. 
spatial perturbations only along the velocity gradient. The linear equations for ^^(k, t) are 
then of the form 



at^„(k, t) + M„^(a, k)^^(k, t) = 



(7) 



with 



M(a,k) 



nik 

cio^ {r]/2e)a^ + Dk^ -{2r]T/e)aik {pT/e)ik 



—Undik —{v /2T)aik 
{pn/p)ik {p/pT)ik 









a 
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where ci = —pUnT/e, D = XT/e, u = rj/ p, = + n/p, and Zn = dz/dn. 

The dispersion relations obtained from det(sl — M) = give five hydrodynamic modes, 
Said, k). If the real parts of one or more modes become negative for some values of k and a 
then solutions to Eqs. (|^) grow in time and the uniform shear flow state is linearly unstable. 
Direct calculation shows this is the case for sufficiently small k for any fixed and finite a. The 
essential mechanism for the instability follows from the linearization of the viscous heating 
term PijdiUj in the equations for 6T which is marginally controlled by the linearization of 
the thermostat. The hydrodynamic modes for finite a and asymptotically small k can be 
calculated explicitly to order fc^. Four modes vanish as — * while one is finite. 



So vk^ , si —> {r]/2e)c? , S2 — > C2k'^ 
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S3 = S4 ^ ick — [di — ){k/ay 



(10) 



where pc^ = [2p{2 + nvn/^) + npn]-, pc^C2 = n{3iypn-2piyn), {pvYdi = p [2e{2v + nvn) - ^p], 
and 2c?2 = (3z/ + 7 — C2). These are positive constants depending only on the density and 
temperature of the reference state. The first three modes are stable whereas the complex 



conjugate pair S3 and S4 are unstable for a? < di/d2. This is a primary observation of our 
work: within the hmitations of the well-estabhshed Navier-Stokes equations (e.g., small k 
and small a) the equations are unstable at all k for a reference state with sufficiently small 
shear rate. Conversely, a similar expansion of the eigenvalues in powers of a at fixed k 
shows the instability at all shear rates for sufficiently small k. These asymptotic results 
are confirmed by an exact evaluation of the eigenvalues from (H), using the Percus-Yevick 



approximation for the virial equation of state [|r^ and using the Enskog kinetic theory to 
specify transport coefficients [Q. The results are shown in Figure ^ indicating the line in 
the k-a plane separating stable (above) from unstable (below) dynamics at three different 
densities (n* = na^ where a is the hard sphere diameter). The instability is qualitatively 
the same for all densities. 

The Navier-Stokes approximation for the heat and momentum fluxes requires k <^ inverse 
mean free path and a <^ inverse mean free time. Since the reference state is generated by the 
shear it is possible that higher order contributions in a might remove the long wavelength 
instability. To address this question, it would be desirable to derive the linear hydrodynamic 
equations from the Boltzmann-Enskog kinetic equation without this limitation to small a. 
However, no solution to this kinetic equation is known, for either the reference uniform shear 
flow state or deviations from it. Consequently, we consider the case of low density for which 
the Boltzmann equation applies. While no solutions to this kinetic equation are known 
either, it is well-established that closely related "kinetic models" provide practical and accu- 
rate representations of solutions to the Boltzmann equation. Here, we chose the non-linear 
BGK kinetic model obtained by replacing the Boltzmann collision operator with an average 
collision frequency times the deviation of the distribution function from a local equilibrium 
distribution. The parameters of this local equilibrium distribution are chosen to enforce the 
exact conservation laws. As a consequence, the Navier-Stokes hydrodynamic equations ob- 
tained from the BGK model are the same as those from the Boltzmann equation - only the 
values of the transport coefficients differ. On this basis, we expect that hydrodynamic equa- 
tions outside the Navier-Stokes limit derived from the BGK model should provide a faithful 

6 



representation of those that would be obtained from the Boltzmann equation. These expec- 
tations for the BGK model have been confirmed in the case of shear flow by both analytical 



and numerical comparisons with the Boltzmann equation [0. It is remarkable that the 
BGK equation can be solved exactly to determine the reference state distribution function 
for uniform shear flow, at arbitrary shear rate |T^. Using this known result, the heat and 



momentum fluxes in the conservation laws (|l|)-® can be determined to leading order in the 
deviations of the hydrodynamic fields from uniform shear flow. The resulting hydrodynamic 
equations linearized about the reference state are again valid only to order k'^ but now there 
is no a priori restriction on the value of the shear rate; the form of Eqs. (|^) are unchanged 
but the matrix elements of Mapia, k) are not restricted to order a^. This new shear rate de- 
pendence leads to qualitative differences from the Navier-Stokes equations (e.g., rheological 
effects such as shear thinning, normal stresses). The single parameter of this kinetic model 
is the average collision rate for the Boltzmann equation and all dependence on the interac- 
tion potential occurs only through the temperature dependence of this parameter. We have 
chosen the simplest case of Maxwell molecules, V{r) ~ r~^, for which it is a constant. All 
transport coefficients of these generalized linear hydrodynamic equations can be calculated 
exactly as functions of the shear rate, and the eigenvalues of M„/3(a, k) can be determined 
just as in the case of the Navier-Stokes equations. A long wavelength instability for any 
value of the shear rate is found again, now including values of a well outside the limitations 
of the Navier-Stokes equation, c.f. Fig.[|. We conclude that the instability observed in (p!OD 
is robust and is not an aberration resulting from the approximations and (|]). 

The extension of the hydrodynamic equations to larger shear rates using the BGK model 
allows comparison with Bird's direct simulation Monte Carlo method [|6|,p!5|. There have been 



significant tests of this method for uniform shear flow |T3|. The method is so accurate and 
efficient that virtually all practical applications of gas kinetic theory far from equilibrium 
now use it. We have used a direct numerical solution to the BGK kinetic equation to 
test the stability analysis without the intermediate step of constructing a hydrodynamic 
description. The solution is constructed as follows. First, the volume is partitioned into 



cells and particles are distributed with positions and velocities according to a specified 
initial distribution. Next, at each finite time step r < mean free time, a streaming and 
collision stage are computed. The particles are moved in straight lines to new positions 
at time t + t. For each particle, the probability of a collision is determined as the (local) 
collision frequency times r. If a collision occurs, the velocity is replaced by a random velocity 
sampled from the local equilibrium distribution. This collision calculation is performed for 
all particles and the whole process is iterated for many time steps. 

To confirm our hydrodynamic analysis of the instability based on the BGK equations 
an initial state for the hydrodynamic variables is chosen with k = 0.1 and a = 0.5 (in 
units of the inverse mean free path and time, respectively). This corresponds to conditions 
for which the linear hydrodynamic equations are unstable. The two unstable modes have 
complex eigenvalues so they lead to an oscillatory time dependence with growing amplitude. 
The hydrodynamic fields are constructed directly from an expansion in the eigenvectors of 
MQ,^(a,k). To illustrate this dynamics we have chosen initial conditions such that 6Uy{0) 
couples only to the unstable modes, while 5f/^(0) couples to the both stable and unstable 
modes. Figure |^ shows a comparison of the results for 6Ux(t) and 5Uy{t) as a function 
of time with those obtained from a Monte Carlo simulation of the BGK kinetic equation 
using the Bird method. The good agreement shows that the instability is not a consequence 
of the assumptions behind the hydrodynamic calculations, and also that these equations 
provide an accurate description of the initial stages of the instability. We have performed 
the simulations for times an order of magnitude longer than that shown. The maxima and 
minima continue to grow, although differences between the theory and simulation become 
more significant. This is expected since the linear hydrodynamic analysis is limited to small 
amplitudes. 

It remains to understand the consequences of this long wavelength instability. We plan 
to extend the low density Monte Carlo simulations to very long times to see if a final 
stationary state is regained. At high densities molecular dynamics simulations appear to be 
stable, except at very large shear rates where a transition to an ordered state occurs It 
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is likely that the long wavelength instability considered here has not been seen due to the 
finite system sizes considered, i.e. k > 27t/L [|T6|. We plan to explore longer wavelengths at 
high densities using both molecular dynamics and an extension of the Bird method to the 



dense fluid Enskog equation fTf^ . 



9 



FIGURES 

FIG. 1. Critical lines for stability determined from Eq. (^) at n* = (solid curve), n* = 0.2 
(dotted cm've), and n* = 0.4 (dashed curve). Also shown are the results from the BGK kinetic 
model for n* = (dash-dot curve). All are in units of the mean free path and mean free time. 

FIG. 2. Time evolution of (a) 8Ux{t) and (b) 5Uy{t) for a = 0.5 and k = 0.1 in units of the 
inverse mean free time and mean free path. 
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